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Abstract 

We study how the spatial distribution of inertial particles evolves with time in a random flow. We 
describe an explosive appearance of caustics and show how they influence an exponential growth 
of clusters due to smooth parts of the flow, leading in particular to an exponential growth of the 
average distance between particles. 
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Random compressible flows generally have regions where contractions accumulate and 
density grows. Infinitesimal elements expand or contract exponentially which can be char- 
acterized by the set of Lyapunov exponents. Since the sum of the exponents is non-positive 
, density tends to a singular multi-fractal set with moments growing exponen- 
tially. Both the evolution and the final state of density in spatially smooth random fiows 

~\ . The fiow of inertial particles is 



have been described recently within some models yl y, 
compressible even when the fiow of ambient fiuid is incompressible j7| so particles participate 
in the fractalization and have some of their concentration moments growing exponentially 
0]. On the other hand, every time there is a negative velocity gradient exceeding the inverse 
viscous response time of particles, faster particles from behind catch slower ones creating 



folds in distribution and caustics 



Such breakdowns of distribution lead to finite-time 
singularities and explosive growth of some density moments. The goal of the present paper 
is to describe the statistical evolution of concentration from a uniform one to a set of clusters 
and voids and, in particular, to describe the role of of folds in this evolution. 

Because of folds, the problem of inertial particles in a fiow is kinetic rather than hydrody- 
namic Analytic approach to a realistic kinetic description does not seem to be feasible 

now. On the other hand, the significant progress of analytic Lagrangian methods makes 
it tempting to use them: to follow, for instance, a couple of close particles and to account 
only for a local velocity gradient. The question is: what can we learn from the Lagrangian 
approach about the statistics of particle concentration? To answer that, one needs a model 
that allows to compare numerical data from kinetics with an analytic Lagrangian solution. 
For that end we consider here the motion of inertial particles in a one-dimensional random 
fiow, which is appropriate for our main goal to describe the role of breakdowns that are 
one-dimensional in any space dimensionality. This model is a subject of much interest from 



different perspectives 



Here we briefly review what is known and derive new results. 



in particular, describe the statistics of the inter-particle distances R. We also carry direct 
numerical simulation of kinetics in this model and flnd the growth rates of the moments of 
concentration n. It is only for smooth flows that one can immediately convert R into n (in 
Id simply taking n = Since the flow of inertial particles has discontinuities, any 

given interval between two chosen particles does not contain the same particles all the time. 
Particles can enter and leave the interval i.e. numerous folds appear in particle distributions 
making nonlocal even the problem of describing single-point density statistics. We show 
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that indeed the growth rates of density moments and inter-particle distances are different. 

Particle coordinate q and velocity V change according to dq/dt = V(q, t) and dV/dt = 
[u(q, t) — V]/r with q(r, 0) = r. Here the viscous (response) time is r = {2 / 9) {po / p) {a^ / u) 
with a particle radius and pQ,p particle and fluid densities respectively. We treat the fluid 
velocity u as a given random function of time and smooth function of space coordinates. Let 
us briefly remind some relevant properties of smooth compressible random flows The 
behavior of an infinitesimal volume is governed by the local matrix of derivatives (called 
strain matrix) taken in the Lagrangian frame Sik = dui/dxk- Considering the distance 
between two fluid particles, R(t, ri — r2) = q(ri, t) — q(r2, t) one finds {R"^) ~ exp(£'mt) with 
being a convex function of m. Density can be expressed as n(t) = dei'^ dRi{t^Y) / dvj 
(provided that the initial distribution is uniform tiq = 1) so that the Lagrangian moments 
(n"™") are related to space-averaged (Eulerian) moments via (n""*) = {n^~^)E ~ exp(rmt) 
(every trajectory comes with the weight n~^). Therefore, Fq = = Fi which correspond 
to conservation of mass and volume (Lagrangian and Eulerian measures) respectively. In 
one-dimensional (Id) smooth flows, F^ = Em- 

In Id, one has for the distance R{t) and velocity difference v{t) between two close inertial 
particles: 

R = v, Tv = sR-v tR + R = sR . (1) 

The substitution R = \1/ exp(— t/2r) turns into Schrodinger equation with a random 
potential (Anderson localization), with space and the localization length replacing time and 
the Lyapunov exponent. 

The quantity a = v/R satisfies the Langevin equation driven by the random noise s{t) 

& = -a"^ - a/r + s/t = -dU/da + s/t . (2) 

Let us describe the probability of finite-time singularity (explosion) cr — — oo which corre- 
sponds to crossing of particle trajectories. Such probability can be written as a path integral 
over trajectories with cr(0) = (Tq? c"(T) = — oo: 

rT 



P{T)= J VaVpVsV{s}exp(^J^ 



& + W-- 

T 



dt' } . (3) 



Here V{s} is the probability functional for s and W = U' = {a"^ + a /t). When T is less than 
the average time between explosions (defined below), P{T) is given by the single trajectory 
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(optimal fluctuation 



12 
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Ql) which maximizes the probability and can be found by a 
saddle-point integration of 

First, consider T which is much less than the correlation time of the air gradient s. Then 
the optimal fluctuation corresponds to s = Sq which does not change during t. In this case 
the integration over the processes s{t) is reduced to the averaging over a single value Sq with 
the measure Ps(so), which is a single-time statistic of velocity gradient s. The saddle-point 
integration over the fields p, a is reduced to solving the equation Q with constant s{t) = sq 
and the boundary conditions a{0) = ao, cr(T) = — oo. Straightforward integration yields the 
following relation: 

dcr 



T 



T 



So 



a — ra^ 



(4) 



r(-l-4sor)-i/2 



TT — 2 arctan 



2(ToT 



V-1 -4sor^ 

which formally gives a relation between the optimal value of Sq and the collapse-time T. It 
is not possible to find the analytic expression for so{T) for a general value of do, however 
the situation greatly simplifies for o"o = -|-oo. In this case the PDF P{T) can be interpreted 
as the distribution of time intervals between consequent collapses. Note, that as long as the 
trajectory starting from ctq = +oo passes through all values of a this distribution is also a 
lower estimate for the P(T) for general (Tq. Substituting do = -|-oo in (jH) one obtains 

~ (5) 

(6) 



or equivalently 

"47 1^ 

In this case the probability of collapse is given by 



So 



P{T) = Psiso) 



dSr 



dT 



27rV 



1 

'47 



(7) 



T3 ' V 4r T2 

One can see from this expression that collapses occur only if there is a finite probabil- 
ity of having sufficiently negative flow gradient, s < — l/4r. In particular for Gaus- 
sian gradients, Ps{x) = (a/vr)-'^/^ exp(— ax^), the short-time asymptotics is as follows: 
P(T) ~ T-3exp(-a7r2rVT4). 

Consider now the case when the correlation time of s is much shorter than T. In this 
case, the noise can be effectively considered as white Gaussian, {s(t)s{0)) = 2DT'^6{t), and 



P(T) = I Vaexp!^-^l\a + W]'dt'^ 



(8) 



For DT^ -C 1, if follows from the saddle point approximation that the probability is given 



m 



14| ) which satisfies 



by the optimal fluctuation (also called "instanton" trajectory 
a = W{a)W'{a) with the boundary conditions a{0) = ctq, cr(T) = — oo. After one integration 
one obtains the following equation: 



a 



-VE + (9) 



where E is an integration constant, characterizing the trajectory. This constant is deter- 
mined by the boundary conditions: 

n da , , 

T= / , 10 

7-00 ^E + W 

The probability of such fluctuation is given by P(T) ~ exp(— A), where 



Unfortunately, the integrals p()|lHl can not be expressed through known special functions, 
so we are able to get analytical results only in some limiting cases. We will consider the case 
(Jo = +00 following the same arguments as in the preceding analysis. First, we consider the 
limit Et^ <^ 1 which as follows from ()10|) corresponds to large times T ~ rlog(l/£'r^) ^ r. 
From the expression (|TT|l we have in the main order: 

p {\W\-W)'da ^ f \W\da ^ J_ 

loo 4D|l^| D QDt^- ^ ^ 

We see, that in the main approximation the action does not depend on the T, which has 
a simple interpretation: the collapses are produced by universal tunneling processes, each 
having a probability exp(— l/6Dr'^) and characteristic time-scale r. In order to find the T 
dependence of the total probability we should study the fluctuations around this instanton 
which would involve some bulky calculations. However, for the intermediate region of 
T <^ r exp(l/6Dr^) one can treat these tunnelings as a Poissonian process and predict the 
linear behavior P{T) ~ T/r exp(— l/6i5r^). This expression is certainly not true in the 
case Dr^ ^ 1 when the action A is not large and the saddle-point approximation is not 
applicable. Another limiting case, which can be studied analytically corresponds to the very 
high "energies" Et ^ 1 where one can neglect the linear cr/r terms in (jlOfllll . so that one 
has 

^ nm' _ r(i/4)^ 31.5 



The crossover between the two regimes happens at T ~ r. To summarize, for the white s{t) 
one gets 

^ ^ f exp{-c/DT^), T <T, 
P{T) { ' " (14) 

I Texp(-l/6L'r3) r <T < t exp{l/6DT^), 

where c = [r(l/4)]V967r2 ^ 31.5. 

Since we consider dilute distribution of particles and neglect their pressure, then a changes 
sign after the explosion as the fast particle overcomes the slow one. That is the flux of 
probability that goes to a ^ — oo returns from a —>■ +00. That allows for the steady-state 
probability density function (PDF) having constant probability F flux equal to the number 
of breakdowns per unit time. Such PDF must have -P(cr) ~ Fa^^ at a ^ ±00. If, as 
is usually the case, the initial P(cr, 0) does not have power tails, they appear at t = +0 
according to P{t,a) oc P(t)cr"^ and d7ITi|) . 

When a —00, i? — 0. To establish the sufficient condition for negative moments of 
the distances to blow-up in a finite time, introduce Ri^k = {c^R^)- Assuming even fc, using 
(j21) and Cauchy inequality Ri^k < ^Jk^ol S^^ ^"^^ ^ ~ majoring inequality 

k{Ztt + Zt/r) > 0. For positive k, it means smooth evolution with Z growing. For negative 
k, this inequality gives Z{t) < Z{ti) + TZt{ti){\ — e"'-*"*^^/^). This means that Z turns into 
zero, and respectively, the negative momenta of the distances {k < 0) will blow up in a finite 
time if at some ti: Z + rZt < or in other terms, rdRQ^k/dt > |/c|i?o,fc- This condition 
is readily satisfied for most random processes s{t), the detailed analysis will be published 
elsewhere. 

In the rest of the paper we approximate the fiow gradient s{t) in the particle reference 
frame by a white noise, which is quantitatively good for heavy particles and give a qualita- 
tively correct description in other cases. In the white case, a variety of analytic results can 
be obtained, some translated from the localization theory and super-symmetric quantum 



mechanics 



3,C 



found explicitly [ij 



and some original that we derive here. The steady-state PDF can be 



da' , (15) 



F 


r ^(^)i 






^exp 


Jexp 


D 


D 



with the fiux F = DdPo/da + {a^ + a/r) Pq ^ {27it)-^ exp[-l / (QDt^)] for Dt^ < 1 (the 
dimensionless Stokes number Dr^ = St measures the inertia of the particle). At St ^ 1, 



average time between breakdowns is much smaller than r 

n 

in this limit. The Lyapunov exponent (a) changes sign at St* ~ 0.827 [ll|: (a) ~ —Dt'^/2 
at 5t < St* and (a) ~ D^^^ at 5t > St^. That means that small particles cluster while 
large ones mix uniformly. 

Note that the Gibbs state exp{—U/D) is non-normalizable in this case. The flux state 
(|T5|) minimizes ent ropy production [l2|. It can be shown that it is indeed the asymptotic 
solution at t — s> cx) 

To describe the joint statistics of a and R we introduce the generating function Zk{cr, t) = 
{6 [a (t) — a] R^{t)), which satisfies the equation 



dt 



, n, d f a 2 

kaZk + - + ^ + ^ir 
oa \T da 



(16) 



Substitution Z^. = \I/((T, t) exp[— t//2D] turns it into the Schrodinger equation in a double 
well, which has been a subject of numerous works related to tunnelling and instantons 
(see e.g. jl^ 21, Following [3, Q we first find (non-normalizable) solutions 

exp(7fct/r — U I D)fk{a) with fk being polynomials and then the conjugated solutions by the 
method of variable constants. For example, there are steady states Zq = Pq and 



Z,{a) 



ar 



) exp 



Via) 
D 



exp 



D J (1 



da' 



a'T] 



In particular, this solution allows one to obtain the mean velocity difference between two 
particles at the distance a: a J aZi{a,t) dt needed, for instance, to calculate the coUision 
rate. The growth rates of the moments of inter-particle distance can be obtained from (|16|1 
or in a straightforward way by writing 



Ri,k = -lRi,k/r -{I- k)Ri+i,k + l{l - l)DR, 



l-2,k , 



(17) 



where the higher moments are expressed only via lower ones. Assuming that for a given k 
all Ri^k oc exp(7fc)f:) we get for 7^ the {k + 1) -st order algebraic equation. For the second 
moment one gets 72 (72 + t^^) (72 + 2r~^) — AD = which gives 72 ~ 2Dt^ for Dr^ <C 1 
and 72 ^ (4L))^/^ for Lir^ > 1. 

For arbitrary k, we find asymptotics. If Dr^k^ -C 1 then 7^ ^ Dr'^kik — l)/2. When 
Dr^k^ ^ 1, the determinant of (|17p is approximately 7^"*"^ — 'j^~'^Dk(k — 1) where 

= YliK^ — i) (X k"^ and 7fc oc {Dk'^Y^^. Let us compare the growth rates of the distance 
moments for the inertial particles with those for smooth compressible short-correlated flow. 
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For the latter, 7^ ~ k{k — 1) while for the former the dependence is parabolic only for 
low-order moments in the low-inertia limit Dr^k"^ <^ 1. High moments correspond to high 
inertia and have 7^ ~ (Dk'^Y^^ even for St <ti 1. Note that conservation requires 70 = 7i = 
for inertial particles as well. 




FIG. 1: Growth rates of distance moments for a smooth flow (broken line) and inertial particles 
for different Stokes numbers (two solid lines). 

Since R is sign-changing for inertial particles, the statistics of |-R| deserves separate study, 
particularly for comparison with the concentration. The equation for the time derivative 
of Rik = {(T^\R\'") differs from jllj) by the extra term 2{a^+^R''+^6{R)) , which IS nonzero 
for I = k. As a result, the growth rates 7^ = Rj^dRik/dt differ remarkably from 7^. The 
most dramatic new effect can be readily appreciated since 7^ are related to the Lyapunov 
exponent via (cr) = {d%/dk)k=o- At high inertia, when St > S't* and (cr) is positive, 
it is thus evident that 71 > as seen from the sketch in Fig. ^ Nonzero growth rate 
of is a remarkable qualitatively new effect with a clear physical meaning: in every 

breakdown, extra particles enter the interval between the two particles that we follow; the 
interval length must grow to ensure conservation of the total number of particles. From this 
interpretation, it is clear that the growth rate must be nonzero at low inertia as well, when it 
must be proportional to the exponentially small rate of explosions: 7^ ~ F ~ exp(— 1/65^). 
Remarkably, one can also establish asymptotically exact pre-exponential factor. Consider 
d\R\/dt = <7\R\ + 25{R)aR^. The growth of (|-R|) must be determined by the last term, 
which accounts for the breakdown processes, since (i?) does not grow. In order to obtain the 
explicit expression for 71 we first analyze the dynamical equation on the stages between the 
breakdowns, which formally coincides for both R and \R\ and then account for breakdowns 
explicitly. For delta-correlated s, we can break the time interval into pieces with independent 
evolution (markovian property). Using this fact and multiplicative nature of (Q) one can 
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derive the following identity 



/\m\_\ ^ / \Rit)\\^/ \v{u)\ \ 
\\Rm/ \\Rm/\HtN)\/l\\\v{h--i)\/ ^ ^ 

Here ti..tAr are the times when the breakdowns happened, v{tk) are absolute values of the 
velocities in these breakdowns. All the averages in this expression correspond to the dynam- 
ics between the breakdowns, for which (Q) can be solved explicitly. In order to study the 
dynamics between breakdowns we introduce the function Z{a, (Tq, t) which is the solution of 
()16|) with k = 1 and the initial condition Z{a, ao, 0) = (5(cr — (Tq). In contrast to the previous 
analysis we consider different boundary conditions for the function Z{a,(7Q,t). Namely we 
will assume that there is no flux at cr = +00, which means that Z{a, (Tq, t) decays exponen- 



tially there. In this case Z{a,ao,t) can be interpreted as Z{a,ao,t) = R{t)/R{0) where the 
averaging is performed only on the trajectories which had no breakdown up to time t and 
which satisfy the following boundary conditions a{t) = a, a{0) = a^. Note, that for such 
trajectories we have \R{t)/R{0)\ = R{t)/R{0). We are able to fix the breakdown moment 
at t by taking the limit cr —00. Analogously the limit ctq — ^ +00 fixes the preceding 
breakdown moment at t = 0. In order to analyze the product (fTHj) we introduce three new 
functions: 

J+iao,t) = -a Z(cT, cro,t)U^-oo (19) 



\R{tk-t)\ 
\Ritk + t)\ 



J-{a,t) = crQ^Z{(T,(To,t)\a^,^oc (20) 



v{tk)\ 

M(tk+l — tk) = CTq V+(cro, tfe+l — tk)\ao~^oo (21) 



1^(4)1 

These function have the following meaning. J+(cro, t)dt is an average ratio of \v{t)/R{0)\ for 
trajectories with the boundary condition a{0) = ctq which had breakdown at the interval 
(t,t + dt) . Analogously J_{(T,t) is a ratio of \R(t)/v{0)\ for the trajectories which emerged 
after breakdown at t = 0. Finally, M{t)dt gives us the average ratio of the velocities v{t)/v{0) 
between the two breakdowns which happened at time and in the interval (t, t + dt). Note 
that the normalization factor in the definition of J^{aQ,t) accounts for the fiow of 
trajectories given by (cr^ + aT~^)Z. With such a normalization one has J+(cro, t) is the 
average ratio of \v(tk)/ R{0)\ of all trajectories which satisfy cr(0) = ctq and had a single 
breakdown at time tk G Q. After introducing these three new functions we are able to 



average the ratio \R{t) / R{0)\ over all trajectories with an arbitrary number of breakdowns 
happened at all possible times. We can write formally 

N N-1 



/\m[\ 

\ 1^(0)1/ 



da Z{a,ao,t) + J2 YldtkJ+{ao,ti)J4a,tk) JJ M{tj+i-tj] 

N=l k=l 7=1 



(22) 



This expression can be simplified by turning to the Laplace transform representation: 

Where the upper index s corresponds to the Laplace transform of a function: 

F' = I dtexp{-st)F{t) 



(23) 



(24) 

Long time asymptotic of both expressions is determined by the most left pole or other 
singularity of the integrated functions. One can easily note all three functions Z*, and 
have the same poles s = E^.. Therefore the long time asymptotic is determined either 
by Eq or the most left solution of the equation = 1. We will show that in our cases 
the asymptotic is indeed determined by the later singularities. First we want to show that 
M° = —1. The Laplace transform of Z obeys the following equation: 



Ddl 



a 



Z„ = 6{a - o-q) 



(25) 



here we have inserted the initial conditions Z{a,t = 0) = 5(cr — ctq) in the r.h.s. explicitly. 
Remarkably this equation may be rewritten in the divergent form after the substition Z = 
(r^^ + cr)^^n. This fact was probably first noted in [2^. New equation acquires the following 
form: 



-a(T 



^ + 



D 



— 1 



Ddl\u=il + ao)5ia-ao) 



(26) 



r + o" 

In order to find 11° we have to set s = in this equation, after which it can be easily 
integrated: 



Z^{a,ao) 



1 r-' + ap 
D T-^ + a 



U{a\ 



J —00 



(27) 



where a' = min(cr, ctq) and U{a) = (r ^ + cr)^ exp(— cr^/2Z}r — a^/3D). Straightforward 
integration of ^I7\i yields the expected result: 



M° = lim 



-1 



(28) 
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where we assumed the hmit a — oo, ctq — > +00. Returning now to the determination of 
long time asymptotic of \R{o, t)\ we conclude that it is given by the solution of the equation 
= 1. In general case the explicit form of can not be found analytically, so in the 
next consideration we will assume the limit Dt^ ^1. In this case we know that the growth 
rate of |-R| is parametrically small, so that the solution s is almost zero. So, we need to 
know what is the behaviour of M'^ near zero. In order to analyze it we turn the equation 
After the substitution H = U^l'^{a)^ we arrive to the quantum-mechanical problem 

il^ = -s^. if = -Da^ + ^!fe^^-l-2a (29) 

Such asymmetric double-well Hamiltonians have been extensively studied in the literature, 
see e.g. Q, 13] where the spectrum of 11 was analyzed in the limit ^ 0. Omitting the 
details we will just note that in the main order the ground state energy of H is negative, 
with absolute value Eq = —E = —{Dt^/2tt) exp(— l/6Dr'^) while the other energy levels are 
positive and are of order unity (are not exponentially small). From the hermiticity of H it 
follows that the general form of will be the following: 



M' = y — ^ (30) 

where Ck = — cr^f/"'^/^(o")f/~^/^(cro)^fc(cr)\E'fc(cro) taken at limits a —oo,ao — > +00. All Ck 
in the main order are proportional to oc exp(— l/6i5r^) and are thus exponentially small, 
while out of all energies Ek only Eq is exponentially small. Therefore in the vicinity of s = 
only the term with = is relevant. Although we could find the cq ab initio we won't 
do that and will instead use the fact that M° = —1, which immediately yields Cq = E. 
Therefore in order to find the expression for the growth rate of |-R| we have to solve the 
algebraic equation E/{s — E) = 1 from which we finally obtain the growth rate exponent of 
\R\. 

= ST = 2Et = (St/n) exp(-l/65t) . (31) 

This final expression shows that indeed the leading singularity in is determined by the 
solution of M'^ = 1. The growth rate 71 is exponentially small because it is determined by 
the rare breakdown events. Let us emphasize that we have established asymptotically exact 
pre-exponential factor in (jHH). 

We now present the results of numerical simulations of the growth of particle separation 
< > in Lagrangian frame and of negative moments of density < n~'^ > in Eulerian 

11 



frame. The method used to obtain the growth rates is the Multicanonical Monte Carlo , 
a technique of adaptive importance samphng which boosts the probabihty of rare events 
that determine large negative moments. The Lagrangian results were obtained solving (Q). 
The results presented in Figs ]2l3l confirm an exponential growth of (|-R|'^). 
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FIG. 2: The moments < \R\^ > for A; = 2, 3, 4 for St = 0.2. Time is normalized by r. 
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FIG. 3: Because of inertia the modulus of particle separation < \R\ > grows {St = 0.2). 



We also observe an exponential growth of the particle separation, < |_R| >. Figure El 
shows a good agreement between the numerics and the theoretical prediction (j3H) up to a 
fairly large St ~ 0.35. 
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FIG. 4: The growth rate 71 vs the Stokes number. The sohd curve represents theoretical prediction. 



For the ID Eulerian simulations the density field is given by the following expression 

n{x,t) = J dxono{xo)6{x(t\xo) — x) (32) 

where uqIxq) is an initial Eulerian density distribution (which we assume uniform) and 
x{t\xo) is a Lagrangian trajectory of a particle. 

This trajectory is obtained from the system of the ODEs (characteristic equations): 



^x{t\xo) = v{t), x(0|xo)=Xo, 
at 



d , , 



v{t) — u{x{t\xQ) 



(33) 
(34) 



where u{x,t) is the Eulerian Gaussian velocity of the turbulent fiow. 

We assume that it is delta-correlated in time and has a spatial correlation length Ic- 

< u{x, t)u{x', t') >= B{x - x')6{t - t'), B{x) = 5oe-^'/'' (35) 

The specific form of the correlation function B is not important. Eulerian field u{x,t) is 
related to the Lagrangian process s(t) (see Eq. (j2))) via s(t) = du{x{t\xo),t)/dx. From 



this it follows that St 



Di 



(r/2) \B" 



{t/11)Bq. Prior to solving system of 



ODEs ()33|) . ()34|) one has to generate ID Eulerian velocity field u{x,t) with the prescribed 
correlation function (j^S)). The algorithm for this is fairly standard (See. e.g. 12^1). First 
we notice that since the field u{x, t) is delta correlated in time its temporal regularization is 
trivial. Introducing discrete temporal step At at each time step, n, we now need to generate 
spatially distributed Gaussain field Un{x) with the correlation property < Un{x)um{x') >= 
B{x — x') 5rnn- In Order to generate the field Un{x) we utilise the Fourier method. Indeed 
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the field Un{x) can be represented as a following Fourier integral 

oo oo 

Unix) = j cos{2-Kkx)[2E{k)Y''^in{k)dk + j sm{27rkx) [2E{k)Y/\,{k)dk (36) 



where ^n{k) and rin{k) are independent real Gaussian processes with the following properties: 



iUk)) = iVnik)) = ^ ^ 

(37) 

< ^n{k)U{k') > = {Vn{k)Vm{k')) = S{k - k') 6mn 

and E{k) is an energy spectrum of the random field m„, it coincides with the Fourier trans- 
form of the correlation function B{x): 



E{k) = j e^''"'''B{x)dx = Bo v^exp [-nHje] (38) 



We then use a discrete version of (jHUj) : 

M 



Unix) ^ ^/ E{0) Ak^i; + Y^ ^j2E{kj)Ak [^^ cos{27Tkj x) + rj^ sm{27rkj x)] (39) 

i=i 

Here we have partitioned the Fourier space into M intervals, so that the wavevectors kj = 
jAk denote the locations of the equispaced grid points. Variables ^" and 77" form a set 
of independent standard Gaussian variables (mean zero and unit variance). Because of 
the nature of the Fourier method the synthetically generated field Un{x) will contain an 
intrinsic spatial period Xp = {Ak)~^. Naturally one wishes to make it much bigger than the 
characteristic scale of the system L. On the other hand one has to ensure that all we have 
enough harmonics in ()39p to sample the peak of the function E{k). These two requirements 
can be met assuming (/cM)^^ < A A; ^ L~^. 

Once we generated the synthetic Eulerian velocity field u{x,t) we use a method of La- 
grangian markers to obtain the Eulerian particle density at each point (using effectively 
formula (jH^ ). We introduce a chain of Nl representative Lagrangian markers connected 
by some fictitious "strings". Each "string" contains a large constant number of uniformly 
distributed real particles. This number is fixed for each string, it does not change during the 
evolution and is determined by the initial density distribution. During the evolution, the 
strings deform according to the Lagrangian dynamics of the initial markers. In particular 
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the occurrence of explosions in Lagrangian frame corresponds to the formation of folds in 
the chain of markers. In order to obtain numerically the local Eulerian particle density at 
a given point we count the number of strings passing through this point and then for each 
string determine the contribution to the density as a ratio Ni/li where Ni is the number of 
particles in the string and /j is the current length of the string. In Fig. El we plot the first 
four negative moments of n. Similarly to Lagrangian moments, Eulerian moments also grow 
exponentially: < n'^^^ >oc exp(rfct). The table compares and Lagrangian 7^ given by 
(fT7|) for St = 0.1 and St = 0.2. We see that Lagrangian breakdowns (Eulerian folds) violate 
Tfc = 7fc that one would have for a smooth flow. We do not have a meaningful parametriza- 
tion for the dependencies of 7^ — F^ on /c and St. It is likely that rare explosions cannot be 
completely disentangled from the exponential evolution. 



k 


7fc 


Ffc 


% - Ffc 


k 


% 


Ffc 


7fc - Ffc 


1 


0.006 






1 


0.028 






2 


0.158 


0.146 


0.012 ±0.003 


2 


0.274 


0.250 


0.025 ±0.002 


3 


0.393 


0.374 


0.019 ±0.005 


3 


0.643 


0.611 


0.032 ±0.005 


4 


0.695 


0.666 


0.029 ± 0.006 


4 


1.098 


1.054 


0.044 ± 0.008 


5 


L054 


1.012 


0.043 ± 0.009 


5 


1.627 


1.564 


0.063 ± 0.009 


6 


1.459 


1.403 


0.056 ±0.010 


6 


2.223 


2.131 


0.098 ±0.012 



TABLE I: The comparison of Eulerian and Lagrangian growth rates for St = 0.1 (left) and St = 0.2 
(right). 

In ID case there is a very simple way of visualizing the dynamics of the caustics. At 
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FIG. 5: The Eulerian moments < n ^ > for St = 0.2. 
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FIG. 6: The evolution of 30 Lagrangian markers with time. The time is normahzed to r and 
increases from a) to d). St = 0.4 

a given time moment, t, one can plot the final displacement of a particle, x{t) versus its 
initial position xq- The Eulerian density distribution can be obtained by projecting the plot 
onto the coordinate axis x{t). At the time moment t = the curve is just a straight line at 
the half the right angle to the axes. During the evolution it will deform, according to the 
Lagrangian dynamics of individual particles (jSni),(ISll) eventually leading to the formation of 
folds illustrating the nonlocal nature of Eulerian density. 

In Figini we plot three stages of evolution of particle distribution. We take A'^^: = 30 ini- 
tially equispaced Lagrangian markers and follow the evolution of the function x{xo) through 
time for a particular realization of the velocity field. We observe that at the initial stage 
(Figinti') the particle displacements are small so that the density distribution is smooth and 
there is one-to-ne correspondence x(t) ^ xq. FiglHb shows the appearance of the first caus- 
tic (a particle overtakes another). At FiglHfc the folds are more pronounced and clearly 
visible. Finally at large times (FigEld) one can evidently observe the effect of the clustering 
of particles. 
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Let us summarize the peculiarities of the evolution of the distribution of inertial particles 
that distinguish them from smooth compressible flows: 1) Infinite moments of density and 
inter-particle distance may appear non-analytically at t = +0; 2) Average distance between 
particles grows exponentially; 3) Moments of density in the Eulerian reference frame grow 
with the rates not reducible to those of distance moments in the Lagrangian frame. The 
work was supported by the Israeli Science Foundation, the EPSRC and the Royal Society. 
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